function [averageBetas maskedBetas averageBetasMatrix stdBetasMatrix] = extractROImask(subjID, ADAPT_t_thresh, FFAcontrast, analysesDir)
%AdaptID_v1_2_Faces: Call the AdaptID_v1_2  driver for face stimuli
%           [params data stimSet] = AdaptID_v1_2_Faces(expInfo)
%
%           AdaptID_v1_2_Faces(expInfo) initializes the driver for the
%   AdaptID_v1_2 experiment. This file should not be called directly.
%   Instead, it should be called from the AdaptID_v1_2() function.
%           The driver reads the experimental parameters, sets up the
%   screen, load the stimuli, and run the experiment. It also collects and
%   saves the relevant data.
%
%
%   Inputs
%   -------
%   expInfo: contains all the experimental information necessary for
%           execution, including the experimental protocol to be ran
%
%   Outputs
%   -------
%   params: Experimental parameters, read from an external file specified
%           on the protocol file
%   data: Data collected during execution, including subject responses,
%           reaction times, and all presentation times
%   stimSet: Information about the stimuli that were presented
%
%   _______________________________________________
%   by Marcelo G Mattar (09/04/2012)
%   mattar@sas.upenn.edu


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    analysesDir = '/Volumes/cluster/jet/mattar/AdaptID/Analyses/';
end

if nargin < 3
    FFAcontrast = 2; %1 for Faces-Objects, 2 for Faces-Scenes
end

if nargin < 2
    ADAPT_t_thresh = 1.65;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BASIC DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Define some basic directories and files
subjectDir = [analysesDir subjID '/'];
FOBSdir = [subjectDir 'FOBSLoc.feat/'];


% Retrieve the analysis directories
FaceDirectories = struct2cell(dir([subjectDir 'Face*.feat']));
FaceDirectories = FaceDirectories(1,:)';
numFaceRuns = length(FaceDirectories);


averageBetas = cell(numFaceRuns+2,8);
averageBetas{1,1} = subjID;

averageBetas{2,1} = 'Run #1';
averageBetas{3,1} = 'Run #2';
averageBetas{4,1} = 'Run #3';
if numFaceRuns == 4
    averageBetas{5,1} = 'Run #4';
end
averageBetas{numFaceRuns+2,1} = 'Average';

averageBetas{1,2} = '1-back';
averageBetas{1,3} = '2-back';
averageBetas{1,4} = '3-back';
averageBetas{1,5} = '4-back';
averageBetas{1,6} = '5-back';
averageBetas{1,7} = '6-back';
averageBetas{1,8} = 'F-test';

averageBetasMatrix = zeros(numFaceRuns,6);
stdBetasMatrix = zeros(numFaceRuns,6);
maskSizes = zeros(numFaceRuns,1);


% Pre-allocate some important variables
feat_cope = cell(numFaceRuns,6);
feat_varcope = cell(numFaceRuns,6);
feat_tstat = cell(numFaceRuns,6);

gfeat_cope = cell(1,6);
gfeat_varcope = cell(1,6);
gfeat_tstat = cell(1,6);

leaveOneOut_cope = cell(numFaceRuns,1);
leaveOneOut_varcope = cell(numFaceRuns,1);
leaveOneOut_tstat = cell(numFaceRuns,1);

maskedBetas = cell(numFaceRuns,6);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD FFA MASK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if exist([FOBSdir 'stats/lFFA' num2str(FFAcontrast) 'mask.nii.gz'],'file')
    lFFAmask = loadniigz([FOBSdir 'stats/lFFA' num2str(FFAcontrast) 'mask.nii.gz']);
else
    possibleContrasts = [1 2];
    lFFAmask = loadniigz([FOBSdir 'stats/lFFA' num2str(possibleContrasts(~(FFAcontrast == [1 2]))) 'mask.nii.gz']);
end

if exist([FOBSdir 'stats/lFFA' num2str(FFAcontrast) 'mask.nii.gz'],'file')
    rFFAmask = loadniigz([FOBSdir 'stats/rFFA' num2str(FFAcontrast) 'mask.nii.gz']);
else
    possibleContrasts = [1 2];
    rFFAmask = loadniigz([FOBSdir 'stats/rFFA' num2str(possibleContrasts(~(FFAcontrast == [1 2]))) 'mask.nii.gz']);
end

FFAmask = lFFAmask + rFFAmask;
FFAmask = logical(FFAmask);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD ADAPT DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract Adaptation copes, varcopes and tstats for each of the Face runs
for runIndx = 1:numFaceRuns
    for adaptIndx = 1:6
        feat_cope{runIndx,adaptIndx} = loadniigz([subjectDir FaceDirectories{runIndx} '/reg_standard/stats/cope' num2str(adaptIndx+1) '.nii.gz']);
        %feat_varcope{runIndx,adaptIndx} = loadniigz([subjectDir FaceDirectories{runIndx} '/reg_standard/stats/varcope' num2str(adaptIndx+1) '.nii.gz']);
        %feat_tstat{runIndx,adaptIndx} = feat_cope{runIndx,adaptIndx} ./ sqrt(feat_varcope{runIndx,adaptIndx});
    end
end

% Extract Adaptation copes, varcopes and tstats for the group level analysis
%{
for adaptIndx = 1:6
    gfeat_cope{1,adaptIndx} = loadniigz([subjectDir 'Face.gfeat/cope' num2str(adaptIndx+1) '.feat/stats/cope1.nii.gz']);
    gfeat_varcope{1,adaptIndx} = loadniigz([subjectDir 'Face.gfeat/cope' num2str(adaptIndx+1) '.feat/stats/varcope1.nii.gz']);
    gfeat_tstat{1,adaptIndx} = loadniigz([subjectDir 'Face.gfeat/cope' num2str(adaptIndx+1) '.feat/stats/tstat1.nii.gz']);
end
%}

% Extract Adaptation copes, varcopes and tstats for the leave one out runs
%{
for leftout = 1:numFaceRuns
    for adaptIndx = 1:6
        leaveOneOut_cope{leftout,adaptIndx} = loadniigz([subjectDir 'Face_LeaveOneOut/Face_LeaveOut' num2str(leftout) '.gfeat/cope' num2str(adaptIndx+1) '.feat/stats/cope1.nii.gz']);
        leaveOneOut_varcope{leftout,adaptIndx} = loadniigz([subjectDir 'Face_LeaveOneOut/Face_LeaveOut' num2str(leftout) '.gfeat/cope' num2str(adaptIndx+1) '.feat/stats/varcope1.nii.gz']);
        leaveOneOut_tstat{leftout,adaptIndx} = loadniigz([subjectDir 'Face_LeaveOneOut/Face_LeaveOut' num2str(leftout) '.gfeat/cope' num2str(adaptIndx+1) '.feat/stats/tstat1.nii.gz']);
    end
end
%}

for leftout = 1:numFaceRuns
    leaveOneOut_tstat{leftout,1} = loadniigz([subjectDir 'Face_LeaveOneOut/Face_LeaveOut' num2str(leftout) '.gfeat/cope2.feat/stats/tstat1.nii.gz']);
    thisMask = and((leaveOneOut_tstat{leftout,1} > ADAPT_t_thresh), FFAmask);
    maskSizes(leftout,1) = sum(thisMask(:));
    for adaptIndx = 1:6
        maskedBetas{leftout,adaptIndx} = feat_cope{leftout,adaptIndx}(thisMask(:));
        averageBetasMatrix(leftout,adaptIndx) = mean(maskedBetas{leftout,adaptIndx});
        stdBetasMatrix(leftout,adaptIndx) = std(maskedBetas{leftout,adaptIndx});
    end
end

%display('end');






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DETERMINE THE VOXELS TO AVERAGE OVER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%{

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BASIC DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract Adaptation copes, varcopes and tstats for each of the Face runs
for runIndx = 1:numFaceRuns
    % Generate filenames
    copeZIPED = [subjectDir FaceDirectories{runIndx} '/stats/cope' num2str(index_of_1_back_adaptation_in_feat) '.nii.gz'];
    varcopeZIPED = [subjectDir FaceDirectories{runIndx} '/stats/varcope' num2str(index_of_1_back_adaptation_in_feat) '.nii.gz'];
    tstatZIPED = [subjectDir FaceDirectories{runIndx} '/stats/tstat' num2str(index_of_1_back_adaptation_in_feat) '.nii.gz'];
    % Gunzip cope, varcope and tstat
    cope_filename = gunzip(copeZIPED, fileparts(which(mfilename)));
    cope_filename = cope_filename{1};
    varcope_filename = gunzip(varcopeZIPED, fileparts(which(mfilename)));
    varcope_filename = varcope_filename{1};
    tstat_filename = gunzip(tstatZIPED, fileparts(which(mfilename)));
    tstat_filename = tstat_filename{1};
    
    % Load data into matlab
    this_cope = load_untouch_nii(cope_filename);
    this_cope = this_cope.img;
    this_varcope = load_untouch_nii(varcope_filename);
    this_varcope = this_varcope.img;
    this_tstat = load_untouch_nii(tstat_filename);
    this_tstat = this_tstat.img;
    
    % Save maps
    cope_feat{runIndx} = this_cope;
    varcope_feat{runIndx} = this_varcope;
    tstat_feat{runIndx} = this_tstat;
    
    % Delete gunziped files
    eval(['delete ' cope_filename]);
    eval(['delete ' varcope_filename]);
    eval(['delete ' tstat_filename]);
    
    % Mask maps
    maskedBetamap = this_cope(:);
    maskedTmap = this_tstat(:);
    maskedBetamap = maskedBetamap(FOBS_FFA_mask(:));
    maskedTmap = maskedTmap(FOBS_FFA_mask(:));
    
    if (length(maskedBetamap)+length(maskedTmap)) ~= 2*FOBS_FFA_maskSize
        error('Something is wrong. Check the code');
    end
    
    % Save important information about this scan
    FOBS_averageBeta(runIndx) = mean(maskedBetamap);
    FOBS_maxBeta(runIndx) = max(maskedBetamap);
    FOBS_averageT(runIndx) = mean(maskedTmap);
    FOBS_maxT(runIndx) = max(maskedTmap);
    FOBS_coefVar(runIndx) = std(maskedBetamap)/mean(maskedBetamap);
    
end

clear runIndx
clear copeZIPED varcopeZIPED tstatZIPED
clear cope_filename varcope_filename tstat_filename
clear this_cope this_varcope this_tstat
clear maskedBetamap maskedTmap



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GENERATE GFEAT MAPS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sumofoneovervarcope = zeros(cubeSize);
sumofcopeovervarcope = zeros(cubeSize);
for thisRun = 1:numFaceRuns
    sumofoneovervarcope = sumofoneovervarcope + 1./varcope_feat{thisRun};
    sumofcopeovervarcope = sumofcopeovervarcope + cope_feat{thisRun}./varcope_feat{thisRun};
end

varcope_gfeat = 1./sumofoneovervarcope;
cope_gfeat = varcope_gfeat .* sumofcopeovervarcope;
cope_gfeat(isnan(cope_gfeat)) = 0;

tstat_gfeat = cope_gfeat ./ sqrt(varcope_gfeat);
tstat_gfeat(isnan(tstat_gfeat)) = 0;

% Save tstat_gfeat in subject's masks directory
if saveBool
    Face_tstat_gfeat_Struct = FOBS_FFA_Struct;
    Face_tstat_gfeat_Struct.img = tstat_gfeat;
    save_untouch_nii(Face_tstat_gfeat_Struct, [subjectDir 'masks/Face_tstat_gfeat.nii']);
end

clear sumofoneovervarcope sumofcopeovervarcope varcope_gfeat cope_gfeat tstat_gfeat Face_tstat_gfeat_Struct


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LEAVE-OUT-OUT APPROACH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for leftout = 1:numFaceRuns
    averageover = 1:numFaceRuns;
    averageover = averageover(averageover ~= leftout);
    
    % Generate the t-stat map using the leave-one-out approach
    sumofoneovervarcope = zeros(cubeSize);
    sumofcopeovervarcope = zeros(cubeSize);
    for thisRun = averageover
        sumofoneovervarcope = sumofoneovervarcope + 1./varcope_feat{thisRun};
        sumofcopeovervarcope = sumofcopeovervarcope + cope_feat{thisRun}./varcope_feat{thisRun};
    end
    
    varcope_gfeat = 1./sumofoneovervarcope;
    cope_gfeat = varcope_gfeat .* sumofcopeovervarcope;
    cope_gfeat(isnan(cope_gfeat)) = 0;
    
    tstat_gfeat = cope_gfeat ./ sqrt(varcope_gfeat);
    tstat_gfeat(isnan(tstat_gfeat)) = 0;
    
    
    % Generate the mask for calculating average effect
    FaceADAPT_mask = tstat_gfeat > ADAPT_t_thresh;
    ADAPT_maskSize(leftout) = sum(FaceADAPT_mask(:));
    
    
    % Save FaceAdapt mask in subject's masks directory
    if saveBool
        FaceAdaptStruct = FOBS_FFA_Struct;
        FaceAdaptStruct.img = FaceADAPT_mask;
        save_untouch_nii(FaceAdaptStruct, [subjectDir 'masks/FaceADAPT_mask_leaveout' num2str(leftout) '.nii']);
    end
    
    % Calculate average and max effects within FaceADAPT mask
    maskedBetamap = cope_feat{leftout}(:);
    maskedTmap = tstat_feat{leftout}(:);
    maskedBetamap = maskedBetamap(FaceADAPT_mask(:));
    maskedTmap = maskedTmap(FaceADAPT_mask(:));
    
    if (length(maskedBetamap)+length(maskedTmap)) ~= 2*ADAPT_maskSize(leftout)
        error('Something is wrong. Check the code');
    end

    % Save important information about this scan
    ADAPT_averageBeta(leftout) = mean(maskedBetamap);
    ADAPT_maxBeta(leftout) = max(maskedBetamap);
    ADAPT_averageT(leftout) = mean(maskedTmap);
    ADAPT_maxT(leftout) = max(maskedTmap);
    ADAPT_coefVar(leftout) = std(maskedBetamap)/mean(maskedBetamap);
    
    
    
    
    
    % Now, calculate the effects on the intersection of the matrices
    % Generate the mask for calculating average effect
    INTERSECTION_mask = logical(FOBS_FFA_mask .* FaceADAPT_mask);
    INTERSECTION_maskSize(leftout) = sum(INTERSECTION_mask(:));
    
    % Calculate average and max effects within mask
    maskedBetamap = cope_feat{leftout}(:);
    maskedTmap = tstat_feat{leftout}(:);
    maskedBetamap = maskedBetamap(INTERSECTION_mask(:));
    maskedTmap = maskedTmap(INTERSECTION_mask(:));
    
    if (length(maskedBetamap)+length(maskedTmap)) ~= 2*INTERSECTION_maskSize(leftout)
        error('Something is wrong. Check the code');
    end
    
    % Save important information about this scan
    if INTERSECTION_maskSize(leftout) ~= 0
        INTERSECTION_averageBeta(leftout) = mean(maskedBetamap);
        INTERSECTION_maxBeta(leftout) = max(maskedBetamap);
        INTERSECTION_averageT(leftout) = mean(maskedTmap);
        INTERSECTION_maxT(leftout) = max(maskedTmap);
        INTERSECTION_coefVar(leftout) = std(maskedBetamap)/mean(maskedBetamap);
    else
        INTERSECTION_averageBeta(leftout) = NaN;
        INTERSECTION_maxBeta(leftout) = NaN;
        INTERSECTION_averageT(leftout) = NaN;
        INTERSECTION_maxT(leftout) = NaN;
        INTERSECTION_coefVar(leftout) = NaN;
    end
    
    
    % Calculate %overlap between masks
    ROI_percent_overlap(leftout) = INTERSECTION_maskSize(leftout) / (ADAPT_maskSize(leftout) + FOBS_FFA_maskSize - INTERSECTION_maskSize(leftout));
    
end

clear leftout thisRun averageover cubeSize
clear sumofoneovervarcope sumofcopeovervarcope cope_gfeat varcope_gfeat tstat_gfeat 
clear maskedBetamap maskedTmap ADAPT_t_thresh
clear FOBS_FFA_mask FaceADAPT_mask INTERSECTION_mask 
clear analysesDir subjectDir FaceDirectories
clear FOBS_FFA_Struct FaceAdaptStruct
clear index_of_1_back_adaptation_in_feat



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SAVE ALL RELEVANT INFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract relevant statistics within FFA mask
allResults{2,2} = FOBS_averageBeta;
allResults{2,3} = FOBS_maxBeta;
allResults{2,4} = FOBS_averageT;
allResults{2,5} = FOBS_maxT;
allResults{2,6} = FOBS_coefVar;


allResults{3,2} = ADAPT_averageBeta;
allResults{3,3} = ADAPT_maxBeta;
allResults{3,4} = ADAPT_averageT;
allResults{3,5} = ADAPT_maxT;
allResults{3,6} = ADAPT_coefVar;


allResults{4,2} = INTERSECTION_averageBeta;
allResults{4,3} = INTERSECTION_maxBeta;
allResults{4,4} = INTERSECTION_averageT;
allResults{4,5} = INTERSECTION_maxT;
allResults{4,6} = INTERSECTION_coefVar;

allResults{5,2} = FOBS_FFA_maskSize;
allResults{5,3} = ADAPT_maskSize;
allResults{5,4} = INTERSECTION_maskSize;
allResults{5,5} = ROI_percent_overlap;





% Create simplified versions of the general output
onlyAverageBetas = cell(5,(1+numFaceRuns+2));
onlyAverageBetas{1,1} = subjID;

onlyAverageBetas{2,1} = 'FOBS mask';
onlyAverageBetas{3,1} = 'Adapt mask';
onlyAverageBetas{4,1} = 'Intersection';
onlyAverageBetas{5,1} = '% Overlap';

for i=1:numFaceRuns
    onlyAverageBetas{1,i+1} = ['Run #' num2str(i)];
    onlyAverageBetas{2,i+1} = allResults{2,2}(i);
    onlyAverageBetas{3,i+1} = allResults{3,2}(i);
    onlyAverageBetas{4,i+1} = allResults{4,2}(i);
    onlyAverageBetas{5,i+1} = [sprintf('%3.2f', allResults{5,5}(i)*100) '%'];
end

onlyAverageBetas{1,numFaceRuns+2} = 'Mean Beta';
onlyAverageBetas{2,numFaceRuns+2} = mean(allResults{2,2});
onlyAverageBetas{3,numFaceRuns+2} = mean(allResults{3,2});
onlyAverageBetas{4,numFaceRuns+2} = mean(allResults{4,2});
onlyAverageBetas{5,numFaceRuns+2} = ' ';

onlyAverageBetas{1,numFaceRuns+3} = 'STD(Beta)';
onlyAverageBetas{2,numFaceRuns+3} = std(allResults{2,2});
onlyAverageBetas{3,numFaceRuns+3} = std(allResults{3,2});
onlyAverageBetas{4,numFaceRuns+3} = std(allResults{4,2});
onlyAverageBetas{5,numFaceRuns+3} = ' ';

if saveBool
    save(['workspace_' subjID '.mat']);
end

clear numFaceRuns subjID i numFaceRuns

%}
